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ABSTRACT 



Three techniques for finding the eigenvalues and eigenfunctions 
are investigated, A typical problem involves a linear homogeneous 
differential equation with variable coefficients of the form 

P(x) y"(x) + P'(x) y'(x) + M(x) . 0 . (1) 

The functions P(x) and M(x) are functions which are positive, or have 
at most isolated zeroes on the fimdamental interval (0,L); uj is a 
parameter. Appropriate end conditions are specified so that the 
problem is self-adjoint. 

The three methods are: tfevrfcon's method, Stodola's method, and 

'tli0 s d.GX'i.vGd. Gzid. 3, coHi'pv.'tGi? 

solution by each method is included in the paper. A second problem 
involving Bessel's equation of order zero is solved using each method 
and a comparison of the eigenvalues and eigenfunctions is made with 
tabulated values. 

The resiflts indicate that Newton's method is to be prefeired 
usually. 
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I, IMPRQ33TJCTION 



Consider the linear homogeneous differential equation with 
variable coefficients of the form 

P(x) y (x) + P'(x) y’(x) + u)^ M(x) y(x) = 0, xe(0,L), (1) 

-| 

where P(x) is of class C , M(x) is continuous, and both functions are 
positive, or have at most isolated zeroes on the fundamental interval 
(0,L), Also given are end conditions of the form 

= 0 

( 2 ) 

^2 7(^2) + b^ = 0 

The function y and the parameter tu are to be determined. 

An equation of the form (1) may arise as a result of solving a 
partial differential equation of the form 

(p(,) (5) 

o t 

by separation of variables with boimdary conditions corresponding to 
(2). The solution of (5) is assumed to be of the form 

Tj(x,t) = X(t) T(t) ( 4 ) 

and standard techniques for separation of variables [8] are employed. 
Equation (I) may also occur as the Euler equation [1] of an isopermetric 
problem of the form 
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« 



( 5 ) 



J P(x) y'(x)^ dx = MIN 
0 

subject to the constraint 



J M(x) y(x)^ dx = 1 . 

0 

This is equivalent to solving the problem 



r P(x) y* (x) dx/ r M(x) y(x)^dx = MIN. 
0 0 



In this case, the Euler equation is 

~ (Py.) - Fy = [P(x) y*(x)] + A{x) y(x) = 0 



where 



and 



= p(x) y(x)2 

= M(x) y(x)^. 



( 6 ) 



( 7 ) 



( 8 ) 



( 9 ) 



( 10 ) 



Calc\ilus of variations treats the minimization or maximization 
of f-unctionals such as integrals. This paper considers three techniques 
for obtaining the eigenvalues and eigenfunctions of Eq, (l). These 
techniques are Newton's method, Stodola's method, and the Rayleigh- 
Ritz method [3, 2, 5 ]* 

Newton's method may be applied to the solution of either the 
isopermetric problem or the equation resulting from separation of 
variables. The principal problem is to obtain a value of the solution 
to the differential equation which vrill satisfy the end conditions 
imposed by Eq. (2). 
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Consider y as a function y(x,tu) on those extremals [l] from 
(^■] > y-j)* Let the derivative of y with respect to to be 

S y (x.tu) 

S(x,tu) . (11) 

A correctional equation 

y(x2,u^) - yCx^) 

“new = “Old - S(x^,..) (IS) 

is derived in Chapter II and the sequence of steps to obtain the 
solution is discussed there. 

Now consider Eq. (l) v/ritten in the self-adjoint form 

- [p(-) y’W] = M(x) y(x) (13) 

subject to the constraints (2). In Stodola’s method, consider (I 3 ) to 
have the form 

L y = u)^ y • ( 14 ) 

where the operator L is defined 

2 

For the eigenvalues cu^ and associated eigenfunction Xi, 

L Xi = cu^ Xi . (16) 

Tahing into account the end conditions, define an inverse operator, 

L ~^ , then 

L"^ Xi = -^ Xi ( 17 ) 

cu. 
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Since differential equation (I5) and the boundary conditions (2) 
are self-adjoint, there is an infinite set of orthogonal eigenfunctions, 
Xi, and an arbitrary fimction, Y^, satisfying the end conditions (2), 
can be expanded in a series of eigenfunctions 



Yq = a^ + a^ + ♦♦» 

Repeated application of the inverse operator (17) "to Yq emphasizes 
the coefficient of X^j and decreases the relative size of the other 
coefficients. The resulting function Y is the corresponding eigen- 
fimction. Further details of Stodola's method are discussed in 
Chapter III. 

The Rayleigh-Ritz method is a procedure for obtaining approximate 
solutions of variational problems. The procedure consists of assuming 
that the desired extremal, y(x) can be approximated by a linear 
combination of n suitably chosen f\mctions 



y(x) §^(x) + C2 §2^ + ' 



(19) 



The are to be determined to effect the desired minimum. Usually 
the $ i(x) are chosen to meet the boundary conditions for any choice 
of C^, An approximation of y(x) may also be made by spline functions 
of the form 

2 r f/ \2* 

y(x) + C2 X + X 



+ 0^ [{x - o.p] + ... + [(x - 



where the C. are to be determined and the o', are in (0,L). The 
1 X 

function ~ “ (^ “ if ^ ^ zero for x < o?^ . 
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The C's are not independent hut must he chosen to satisfy the end 
conditions. In both the polynomial and spline fmction approximations, 
determination of the is accomplished hy solving an algebraic 
eigenvalue problem [ 4 ]. 

In the following chapters, Newton’s, Stodola’s and the Rayleigh- 
Ritz mentod are discussed in more detail end computer solutions of 
(1) by the three methods are presented. 
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II. KEV/TON«S I-IETHOD 



The problem of minimizing the numerator of Eq. (l.7) with the 
denominator constrained to eqvial one is again being considered. The 
problem may be expressed in the form 

L o ^ n 

I = J P(x) y’(x)^ dx = J f dx = MIN (1) 

0 0 



and 



L L . 

J = J M(x) y(x) dx = J f dx = 1 
0 0 



( 2 ) 



Assume the desired solution y(x) = y (x) has been found, 

y(x^_) + = 0 

&2 + 132 ° 



Then 

(5) 



and 



J f\x,y,y’) dx = 1 , y = y’^ , y' = y^ (4) 

0 



A curve satisfying conditions (5) and ( 4 ) is called admissable. 



A. DERIVATION OP NECESSARY COl'lDITION 

The first necessary condition is derived for the case where 
b^ = 0 = b^ . Consider replacing y (x) by y (x) + eT| (x), where 
Tl(x) is an arbitrary function vanishing at 0 and L and having 
piecewise continuous first derivatives. It is necessary that 




( 5 ) 
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or 



S = J ^0^^^ dx = 0 
■■ 0 



( 6 ) 



for all T] such that 



M 

de 



= J “ J dx = 0 



0 y '^0 y 



(7) 



or 



g = ; M^(x) Tl’(x) dx = 0 



In Eq, (6) and (s), M^(x) is defined as 



X 



= (fh - ; fi . 



* 



Let 



~y’ '^0 ^ y = y » y’ = y*‘ 



g.(x) = M.(x) - M- = M.(x) - r M.(x) dx/L 
1 1 ^av ^ Jq ^ ' 



( 8 ) 



(9) 



( 10 ) 



The above necessary condition may be reexpressed. It is necessary 
that 



S = J SqCx) 11* (x) dx = 0 
0 



for all T| such that 

T, 



le " Iq ^1^^^ dx = 0 



( 11 ) 



( 12 ) 



Since 



Jq ^iav = ^iav Jq = ^^av ° 



Now choose 



1]'(x) = Sq(x) + c g^(x) 



(14) 
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and choose c so that T] (x) is admissahle; 



de 



= J g^(go + = 0 



0 



If (x) 0, Eq. (l6) can he solved for c giving 



2 

C = - J g(^(x) g (x) dx / J g (x) dx 
0 0 






Now consider 
L 

r 

0 

then 



dl dJ 

de de 



= ° g-,(x))(gQ(x) + c g^(x)) dx 

= J + c g.(x)j dx ^ 0 , 



Then 



unless 



f >0 

de 



Sq(x) + c g^(x) 



0 . 



( 15 ) 



( 16 ) 



(17) 



(18) 



(19) 



( 20 ) 



Hence this is a necessary condition, known as the First Necessary 
Condition, or the integrated form of the Euler equation. Writing 
Eq, (20) as 

f°, + c f^, - J (f° + c f h dx = (Mq(x) + c M^(x)) (21) 



av 



0 1 

it is seen that f and f enter only on the combination 



■o ^0 , ^1 

F = f + c f . 



( 22 ) 
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The inxegral form of Exiler's eqiiation [l] 



X 



F , - r F dx = M = CONSTAM? 
y* 'Jq y av 

is obtained from (20) by setting 



M = + c M. 

^av ^ av 



Now let 



c = - cu 



then the Euler Equation 



(23) 



( 24 ) 



( 25 ) 



Fy,y, y"(x) + Fy,^ y'(x) + Fy^ y(x) " = 0 (26) 

becomes 

P(x) y"(x) + p'(x) y’(x) = - M(x) y(x) . ,(2?) 

It may be noted that this differential equation is self-adjoint. 

B. COilPUTATIONAL ROUTINE 

It will generally be necessary to generate the cur\’'e by for\\rard 
nvimerical integration. There are two parameters 

c^ = c = - u)^ (28) 

and 

C2 = y*(o) (29) 

which are needed to determine a curve which is an approximation to 
the solution; an estimate of y 2 (o) is needed in order to start 
integration. 
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In order to apply Newton's method it is necessary to evaluate 
derivatives with respect to these parameters. Let 



u = 



and 



V = 



9 



9 y(x) 

9 C2 



( 50 ) 



(31) 



Then differentiating Eq, (2?) with respect to c^ gives 
P(x) u + P (x) u* - c^ M(x) u = M(x) y(x) 



(32) 



Denote the left side of ($2) as Lu, defining the second-order 
linear differential operator L. Then the equations for u may be 
written 



L u = M(x) y(x) 



\yy j 



and 



u(0) = 0 



u (0) = 0 



The corresponding equations for v are 



(34) 



ft ! 



P(x) V + P (x) V - c^ M(x) V = 0 



(35) 



and 



v(0) = 0 
v'(0) = 1 



(36) 
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or 



L V = 0 



v(0) = 0 
v’(0) = 1 



(57) 



The correctional equation for J is obtained from the approximate 
relations 



by J is meant the computed value. The correctional equation for y 
is obtained from the linear approximation 



where y(L) is the computed value. 

However this problem has a peculiarity that makes it somewhat 
simpler. Because both integrals are qtiadr-atic and homogeneous in the 
pair y(x) and y (x) , the principal problem is to find c^ so that 
y(L) = 0. The constant C 2 is then obtained by substituting the 
solution y into J : 




(58) 



= - J ; 



A y(L) = u(L) a c^ + v(L) a 



(59) 



= - y(L) 



C2 = 1/ 



(40) 



Consider then y as a function of a single parameter 



y = y(x, c^) 



(41) 
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The correctional eauatior. reduces to 



A y = u(L) A c^ (42) 

= y(L)A^(L) 

or 

A U)^ = -y(L)/u(L) (45) 

2 

In the computational routine, an initial estimate of aj is made. 

Then y(x) and u are obtained by integrating Eq, (35) ♦ The correction 

2 2 
for tt) is then obtained from Eq. (43)* If a good guess for to is 

made, this method converges rapidly. If a normalized value of y(x) is 

desired then a new solution may be obtained from 

L ” 

^^^^new " (44) 

To tell if the lowest eigenfimction has been obtained, a check 
is m.ade to see if the eigenfunction obtained has a zero in (0,L), 

If it has no zero in (0,L), the.lovrest eigenfunction has been 
obtained. 
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III. STODOLA'S METHOD 



Consider Eq, (1-1) 

^ (p(x) y’(x)) = -u)^ M(x) y(x), (l) 

subject to the constraints 

y(0) + b^ y’(0) = 0 

( 2 ) 

y(L) + b^ y'(L) = 0 

Stodola's method of finding the lowest eigenvalue and eigenfunction 

will now be developed. The system of Eqs. (l) and (2) is self-adjoint. 

2 

Tr; "this <^-3.0 8 r.T*c: iniTi.ni*tG niinibGi' of* Gi.ir8nv3.1vi^s 

i = 1, 2, ... and associated eigenfunctions X^. 

Consider Eq. (l) as having the foirni 

L y(x) = (u^ y(x) . (3) 

in which the linear differential operator L is defined by 

jr'w) (4) 

and consider any particular eigenfunction Xi. Then 

L Xi = m? Xi 
1 

Xi = iu| Xi (5) 



l" XI = »=" 
1 



Xi 
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The inverse operator, 
the end conditions (2), 



_1 

L , is now defined, 
Consider 



taking into account 



^ (p(x) Xi') = - toj M(x) Xi 



( 6 ) 



and integrate (6) to obtain 

P(x) X^' = “ ‘“i (j Xi dx - c) ; 



(7) 



the constant c is to be determined later. The eigenfunction, X^, 
taay be obtained by dividing Eq, (7) by P(x) and integrating again 



X. = - COT 
1 XL 



J. 



J M(x) Xi dx 



0 



'P(^“ 



dx - ! 



P(^J 



dx 






( 8 ) 



The constants of integration, c and d, may be obtained from the 
constraints in Sq. (2), and Eqs. (7) and (S), Solving these 



equations gives 



L L 

^ r M(x) X. dx b„ r M(x) X. dx 

r 0 ^ ^ *^0 ^ 

0 p(l) 

f dx ^2 ^2 \ 'I 

L^2 Jq p( 3^ -*■ P(i7 ■ i^p(oy J 



- (9) 



and 

d = - b^ d/a^ P(0) (10) 

if a^ ^ 0. If a^ = 0 or P(0) = 0 then c = 0 and a solution for d 
is obtained. 

The set of operations on the right side in Eqs. (7) and (S) define 
the inverse operator, L~ , with c and d given by Eqs. (9) and (IO). 
These operations give 
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f 



( 11 ) 



or 



2 -1 

X. = CUT L X. 
11 1 



L"^ X. = x./uj^ 

1 l' 1 



( 12 ) 



Stodola’s method makes use of the fact that an arbitrary piecewise- 
continuous function which satisfies the end conditions (2) can be 
expanded in a series on the eigenfunctions. The first estimate Yo 
of the first eigenfunction may be chosen as follov/s. Let Yo be the 
first estimate of X.: choose it so that it has no zeroes inside the 

X 

interval (0,L). Then Yo may be considered to be the follov/ing form. 



Yo = a^ X^ + a^ X^ + ... (l3) 

where a^ is not equal to zero. For convenience Yo is normalized; the 
norm of Yo, | |Yo| |, is defined as 



MAX |Yo(x) I = I 1Yo| I . 
0 < X < L 



(14) 



and I Iyo I I is set equal to one. Wow consider 



-1 ^2 
L Yo = X + — X^ + . . . 

tu , U) 

1 2 



(15) 



and 



L To - 2m ^1 ■*' 2m ^2 ' 

“2 



oj 2m w ^m 

= ;;;2i ""i + ^2 "^2 ^3 + 



2m 



(16) 
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From Eq. (l6) it caii be seen that the relative size of the components 

/^1 \2 /^i \ 2 

^2’ "* decreasing by a factor » \air~/ ’ *" respectively. 

2 5 

After a few applications of the operator, L~ , to Eq. (l5)» "the 
leading term val]. dominate. It is convenient to normalize after 
each iteration; let 



Z 

m 





1 



I 

where Ym is the m th approximation to and 



( 17 ) 



/ MZmIi ' ( 18 ) 

When this is compared with Eq# (12), where on the left hand side 

corresponds to Y . and X. on the right hand side corresponds to Y • 
m-1 1 ^ m’ 

2 

it is seen that (i)^ is approximately 

= 1/Ii2„,li • (19) 

The iteration is stopped when 

liy^) - ^ ® ’ (-°) 

where e is some specified small number. The resulting function 
is an approximation to . 

In the computational routine an initial approximation, Y^ to 
X^ is made which satisfies the constraints (2) and Y^ is normalized. 

The next approximation to X^^ is obtained by forvfard integration of Eqs. 
(7) and (s). The contrants of integration, c and d are determined 
by Eqs. (9) and (IO). The value of is approximately 1/||Z^|1. 

The new approximation, is normalized, Eq. (I8), th procedure is 
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repeated. When the norm of the difference between successive iterates, 

1 ’ ^m’ less than some prescribed value the routine is stopped. 

2 

The value of u) from Eq. (19) is taken as the lowest eigenvalue and 
is the corresponding eigenfunction. The numerical solution of 
Bessel’s equation of order zero gave a good approximation to the 
eigenvalue but a smaller value of the eigenfunction. 
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IV. RAYLEIGH-RITZ MEl’HOD 



The Eayleigh-Ritz method has long been used to obtain an approximat 
solution to eigenvalue and eigenfunction problems. This is still 
a useful method if a high speed computer is not available. 

Assume that m fiuictions (x) are given which satisfy the end 
conditions. These are often chosen as polynomials or trigonometric 
functions. The solution then is approximated by a linear combination 
of these 

m 

y = S c. §.(x) , (1) 

i=1 

This is substituted into 

p ^ I p ^ p 

0 )^ IS r P(x) y (x) dx / r M(x) y(x) dx = MIN (2) 

0 0 

and the c's are chosen to minimize (2). This is equivalent to the 
minimization of a quadratic form 



v/here 



m m 



S 

i=1 



1 

i=1 



a. . 






c 



i 



c . 

a 



> 




( 5 ) 



subject to a constraint that 
m m 



2 E 

i=1 j=1 




c . c . 



1 J 



= 1 , 



( 4 ) 
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where 



J ^i(^) ^ • 



Eq-uation (2) may then he vjritten 



,2 . , ,.n^ 

^ B C 



where 



A — 9 i — 1,2,###, m , j — 1,2,###,m, 



B — , i — 1,2,###, m , j — 1,2,###,m, 



and 



C — 



c } 



( 5 ) 

( 6 ) 

( 7 ) 

( 8 ) 

( 9 ) 



The problem may he reduced as follows. Let the eigenvalues of 
2 

B he , i = and the associated normalized eigenvectors 

u . . Then 

X 

T^ B T = I (10) 

where 

T = (u^, u^, u^) DIAG (lA^, IA 2 , ♦.*, (11) 



and the transformation 



where 



C = T D, 



T 

D = (8.^ > ^ 2 * ' * ' > dm) 



( 12 ) 



(15) 
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gives 



— JJi — ,*jp ^ 

D D = D T B T D 

-jp 

= C B C 
= 1 

The desired minimum value of 

— — “{P “{T ^ 

C A C = B T A T B = ICH 



( 14 ) 



( 15 ) 



subject to ( 14 ) is the minimum eigenvalue of 

E = T^ A T (16) 

and the corresponding vector B is the associated normalized eigen- 
vector of E. Since the elements of T and B are known, C may be 
evaluated from (4. 12), 

The minimization technique using the spline furiction approximation 
y(x) + C 2 X + [(x - [(x - 

( 17 ) 

is similar to the polynomial approximation method. In this case 
however, the functions ^^(x) do not satisfy the end conditions. 

Hence the constants above are not independenb but one of them is 
dependent on the others for the end conditions to be satisfied. It 
is necessary to have an initial subprogram to eliminate this constant. 
Then the problem is reduced to the same form as before 

m-1 m-1 

X X a. . c. c . = ICQ-I 
i=1 . , xj X 0 

^ ' 0=1 
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and 



m-l m-1 
E E 
i=1 j=1 



b. . 



c . c . = 1 

1 a 



( 19 ) 



and solved as before. 
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V. APPLICATION OF RAYLEIGH-RITZ , 

STODOLA AI'IL Ng.^ON I'lETHOLS 

In Chapters II, III and IV, three methods were developed for 
finding the lowest eigenvalue and eigenf miction for Eq, (l.l). These 
methods were used to obtain numerical solutions of the following 
two problems: Bessel’s equation of order zero 



X y (x) + y*(x) + 10^ X y(x) = 0 , xe(0,2.405), (l) 

subject to the constraints 

y(0) = 1 , y(2.405) = 0 , y'(0) = 0 ; (2) 

and an equation of the form 

(1+x) y (x) + y’(x) + y(x) = 0 , xe(0,l) , (5) 

subject to the constraints 

y(o) = y(i) =0 ‘ (4) 



In this chapter the computation and numerical results are discussed. 

In the Eayleigh-Ritz polynomial approximation, an approximation 
for the extremal, y(x), satisfying Eq. (3) was made by choosing 
suitable functions (x) and ^ 2 ^^) satisfy the constraints in 
Eq* (4)* The approximation for y(x) is 

y(x) = c^ §^(x) + C2 ^5) 

2 

= C^ X (l-x) + c^ X (l-x) 
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The values of , c^, and the lowest eigenvalues were obtained by 
the techniques described in Chapter II. Computer Output 1. shows 
the values of c^ and c^j the lowest eigenvalue and the values of 
y(x) at tenths of an interval. The numerical solution of Eq. (j) 
using the spline function approximation 

y(x) « c^ X - 1 - C 2 ^ [(x-l/3)“] + c^ [(x- 2 / 3 )^] (6) 



is shown in Computer Output 2. and Computer Program 2. This solution 
followed a format similar to the polynomial approximation solution. 
However, a transformation was required to eliminate the dependent 
coefficient. In this case the coefficient c^ was taken as dependent 
on c^ , C 2 and Cy The req'uired mairix products, transpositions and 
eigenvalv.e solutions v/ere evaliiated by using subroutine GMPEI), GOTHA, 
and EIGEN respectively, [ 4 ]* 

A numerical soltition of Eq. ( 3 ) by Stodola's method \ising as a 
first estimate 



y(x) = 4x (l-x) 



'(7) 



to satisfy the constraints in Eq. ( 4 ) yielded the values in Computer 
Output 3- The solution of Eq. (3) by Newton’s method with an initial 
estimate of 

= 2.0 ( 8 ) 



is shown in Computer Output 4» In these two methods the fundamental 
interval (0,1) was sub-divided into one hundred equal sub-intervals 
to perform the second order Riinga-Kutta integration routine. If the 
error between tlie computed value and the boundary- value at X equal 
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one was less than one ten thousandth, the program was terminated. 

The programs for Stodola's and Nevrfcon's method are shovm in Computer 
Program 3. and Computer Program 4* 

The lowest eigenvalues obtained by the above methods are reasonab].y 
close in value. The maximum difference in values was between Newton’s 
method 



10 = 3.78634 

and Rayleigh-Ritz polynomial approximation, 



(9) 



to = 3.79813 . (1.0) 

Newton's method converged in three iterations and Stodola's method 
required seven iterations. 

The RayD.eigh— Hit” , Stodola and Newton methods were ijsed. to f.ind a 
nvtmerical solution to Bessel's equation of order zero, Eq. (l). A 
tabulated solution, [7]» of Eq. (l) gives the smallest eigenvalue 
as one. The result obtained by the Rayleigh-Ritz method was 

to = 1.01435 (11) 



when the function v/as approximated by quadratic splines 






y(x) = c,, + c^ + c^ |^(x - 0.8)^J + c^ (x - 1.6) 



( 12 ) 



The smallest eigenvalue of Eq. (l) obtained by using Stodola's 
method with an initial estimate of 



y(x) := (d-x)/2.405) 



s2 



(13) 
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uu = 1.00000 



(14) 



is 



In Newton's method an initial estimate of 

io = 1.4 (15) 

yielded a value 

oj = 1.0000 (16) 

for the lowest eigenvalue. 

A comp.arison of the eigenfimction, y(x) , obtained by the three 
methods was made with tabulated values for y(x) of Eq. (l). The 
values obtained by Nevrton's method at tenths of an inteirval v:ere 
w5th.in five ten-thousandths, by Rayieigh-Ritz were within five 
thousandths and by Stodola's method v;ere within tv/o one~hrndr‘edths. 

The best approximation of the lovrest eigenfunction was obtained by 
Newton's method. 

2 

The solution of u) obtained by the Rayleigh-Rits method is 

greater than the theoretical value as given in Reference 5« The 

numerical solution of Bessel's equation of order zero supported this 

2 

conclusion. The value obtained for m was approximately one-thousandth 
larger than the theoretical value. 

The results of the three methods are shown in Computer Output 

5., 6,, and 7. and the programs are shown in Computer Programs 5*> 

6. , and 7. 
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VI. CONCLUSION 



The three methods each yield satisfactory values for the eigen- 
values and the eigenvectors. Generally, Ne\'H:on’s method seems to be 
most satisfactory. It is for the most part straightforward to program. 
The only significant errors that are not apparent are there due to the 
integration routine. It is very flexible, being applicable to many 
different types of problems. It converges rapidly if a good initial 
estimate is made, and computing times seem to be generally small. 

Proof of the exact convergence of Ne\'rton’s method may be found in 
Reference 3* 

Stodola’s method also is quite satisfactory. There seems to be 



no questioii of convergence. The initis.1 gi^cs 



n VC mi 



■1 +Q 



and it will still converge. More iterations than for Newton’s 
method seems to be required to get the same accuracy, and there is no 
simple way to get an estimate of the error. 

The Rayleigh-Ritz method seems to have little to recommend it 
if a good computer is available. It is more difficult to program. 

It has one advantage that no iteration is required. Use of spline 
functions increases the tedium of programming so that they are not 
v/orthv/hile . 

The methods of Stodola and Rayleigh-Ritz were used before the 
advent of large scale computers , though application v/as rather limited 
compared to problems that can be treated nov/. 
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In case larger eigenvalues and associated eigenf;mctions are 
needed, ITevrton’s method yields these directly if the initial estimate 
of u) is changed^ The number of zeros of the eigenfunction obtained 
show which eigenvalue and eigenfunction ha.s been obtained. The 
eigenfunctions are automatically orthogonal with a weight factor 
M(x) since the system is self-adjoint. 

If Stodola’s method is used to obtain higher eigenvalues and 
eigenfunctions it is necessary to adjoin a, further condition that 
each new eigenfunction be orthogonal to all of the proceeding. This 
increases the tedium of programming and decreases the accuracy, 
since preceding values are only approximate. 

In the Rayleigh-Ritz method, an estimate of other eigenvalues 
and eigenfunctions is obtained from the other eigenvectors of the 
matrix E of chapter JY, The eigenvalues are too high unless the 
eigenfunction can be expressed as a finite sum of the functions 
$^(x). 
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COMPUTER OUTPUT 1. 



(Polynoiaial) 



THE SMALLEST EIGENVALUE 3o 798129 

THE VALUES OF Cl AND C? 

Cl= 7oC142C'l 

C2= -3ol385?0 

THE VALUES OF Y(X) AT OME-TENTH INTERVAL 



X 


Y 


OoC 


Oo<^ 


Oc IC'CCCf' 


Oc 6.' 3'' ■='9 


Oo2Forro 


1=C21852 


O,3C0CC’C 


1 2T52 72 


O, ■ 


1 ■ 3821 ?“ 


Oo 


Ic 3 612 56 


^0 599999 


lo 231A82 


Oo 699999 


1,011638 


'^,799999 


0- 720556 


9o 899999 


0 377C68 


Co 999999 


0 000'" ^ 
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COMPUTER OUTPUT 2. (Spline) 



THE SMALLEST EIGENVALUE 3o 790108 



VALUES OF 


C1,C2,C3,AND C4 


Cl = 


6o358921 


C2 = 


-7o 197122 


C3 = 


0o565650 


C4= 


5o281205 


VALUES OF 


Y(X) AT ONE-TENTH INTERVAL 


X 


Y 


OoO 


OoO 


Oo ICCOOO 


Oo 563920 


Oo 200000 


Oo 983899 


0o3C0CC0 


lo 259933 


OoAOOOCO 


lo394543 


OoSCOCOO 


lo 395892 


Do 599999 


lo 264613 


Co 699999 


lo 006572 


Oo 799999 


Oo 698055 


Oo 89999 9 


Oo 362533 


Oo 999999 


Oo 000003 
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COMPUTER OUTPUT 3. (Stodola) 



NI= 7oOMEGA = 3078-^142 

THE MAXIMUM DTFFEREMCE BETWEEN YB AND YBN IS OoCPCf-47 
PE,KE,AND ROOT (PE/KE) ARE 7.172C87 Co498484 3o7Q3123 

X YBN 



Oo'' 


OoC 


Oo lOOCC^'- 


387621 


2C0OC'"' 


?o 6942 88 


^0 3CCCor 


Oo 897866 


Oo4cccrc 


<^0992039 


''■c sroccn 


Oo 982383 


Oo 6CDOrO 


Oo 383094 


r, 70000^ 


o-o 713473 


'“V soooro 


Oo 495 322 


Oo9orcco 


Oo2 5f 609 


loOOOOCO 


1 

0 

0 
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COMPUTER OUTPUT k. (Newton) 



THIS IS 


PATH NOo 3o 0^1EGA= 3,78634 


X 


Y 


C 

o 

C 


OoO 


OolOCOO-' 


0o093143 


Co2corot 


C*ol66838 


3O0C0T 


To?15768 


Cc 4CCCCC 


Oo 238397 


CoSCCOOr 


0,236090 


^ ,6C0C0C 


9,212242 




'“o 171498 


'-0 BO or no 


Oo 119087 


Oo 9^0000 


Oon6''29n 


loTOCcro 


9,00*^0 81 
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COMPUTER OUTPUT 5- (Spline) 



THE SMALLEST EIGENVALUE loOOOlCB 



VALUES OF 


Ci,C2,C3,AND C4 


Cl = 


-lol32401 


C2 = 


0o272184 


C3 = 


-O 0 I 25127 


CA= 


-Oo 184540 


VALUES OF 


Y(X) AT CNE-TENTH INTERVAL 


X 


Y 




-lo 132401 


Oo 2 AO 500 


-lo 116657 


OoARlOOO 


-lo069427 


0 0 7 2 1 5 0 i’ 


-Co990712 


Oo 962000 


-Oo 883794 


lo2C2499 


-Oo 759094 


loAABOOO 


-0o617383 


lo 683499 


-Co 459947 


lo 924000 


-Oo 3U2293 


2o 164499 


-Oo 148986 


2 0 405000 


-O 0 OOOOC 9 
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COMPUTER OUTRIT 6. (Stcdola) 



NI= 7o CMEGA= loCOCOOC 

THE MAXIMUM DIFFERENCE BETWEEN YB AND YBN IS OoOOOC21 
PE,KE,ANO ROOT (PE/KE) ARE Co755341 0o761183 Co996155 



X 


YBN 


O 

o 


Ic 000 COO 



0o2405CC 0o974024 



Co 481000 


0o931927 


Oo 721 50C 


Oo 863787 


Oo 962000 


Co772530 


lo 2C2499 


Oo 662059 


lo 44 30 CO 


Oo 537C61 


lo 683499 


0o402789 


lo 9240( 0 


Co 2 64 80 8 


2o 164499 


Oo 128741 



2o4C5C00 OoO 
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COMPUTER OUTPUT 7. (Neivton) 



THIS IS PATH NO. 4 



X 

0.0 

0.240500 

0.481000 

0.721500 

0.962000 

1.202499 

1.443000 

1.683499 

1.924000 

2.164499 

2.405000 



Y 

1.000000 

0.985594 

0.942999 

0.874051 

0.781713 

0.669932 

0.543451 

0.407585 

0.267966 

0.130282 

0 . OOOOiO 



OMEGA = 



0.99998 



COMPUTER PROGRAM lo 



P.AYLEIGH-RITZ METHOD APPLIED TO EQN (lol) 
d + X)*Y' • ( X)+Y’ (X) + ( X)=C, 



DIMENSION dh(2,:..^i ) ,PHP(2, i:i) »A(2,2I ,B(2,2),Pa"’l) , 

1 FA( r 1 ) , FBM ■ 1 d BL( 2 ,2 ) » AM3(2 , 2 ) , C ( 2 ) , AX ( 3 ) t BX ( ? ) , 
2RLB( 2) , BXD( 2) » T( 2, 2 I ,D( ? ) ,TT( 2,2) ,ALPHA(2 ,2) , AP(2 ,2) , 
3S(2,2) .RLBD(2) , Y (ICl ) , Z ( 10 1 ) 

REAL LEV 
X=C-o' 

DX=D-'' 1 
DO 1 I=l,i'l 

C DEFINITION 0^ ELEMENTS 0?= MATRIX 3 
P( 1)^=1 C+X 
PH( 1 , T )=x*( ur-x) 

PH (2, I ) = (X**2) -M loD-X) 

C DEFINITION OF ELEMENTS OF MATRIX A 
PHP( 1, I ) = l, '-2 2^-x 
PHP(2, I )=2o'‘=^'X-?-.C-(X-<=-X2) 

X=X+D„- 1 

1 CONTINUE 

C COMPUTE A(I,K) AND 8(I,K) 

DO 2 1=1,2 
DO 3 K=l,2 
FAINT =•' 

FBINT 

FA ( 1 ) = P { 1 ) -^PMP (1,1) "■•PH" ( K , 1 ) 

F B ( i ) = P H ( I 4 1 ) '■ PH ( K. , 1 ) 

DO J = 2,ri 

FA( J )=P ( J ) -PHP( I , J )*PHP( K, J ) 

FB( J)=PH(I , J)>.'PH(K, J) 

FAINT = FA!NT -+-(FA(J)+ F A( J-1 ) ) ''DX/P 
FBINT = FBINT+ ( F B ( J ) + FB ( J- 1 ) )-' DX/2 
4 CONTINUE 

A( !,K)=FAINT 
B( I,K)=F8INT 
3 CONTI NUF 

2 CONTINUE 

C INITIALIZE BX FOR THE INPUT TO EIGEN ROUTINE 
NC=1 

DO ?''l J = l,2 
DO 2.' 2 1=1,2 
BX(NC ) = B( I , J ) 

NC=NC + ] 

IF (loGEoJ) GO TO 2-'-l 
202 CONTINUE 
2'"1 CON'^IMUE 

CALL EIGEN (BX,S,2,'“) 

C RLB(I) IS 1/SQRl OF EVAL 0^ B 
KK = 1 

DO 8 11=1,2 
BXD( I I ) = BX(KK) 

KK=KK + I 1+1 

8 continue 

DO 84 I 1=1,2 
DO 9 JJ=1,2 
B ( 1 1 , J J ) = ' ^ 

9 CONTIN'T: 

84 CONTINUE 

DO 6 1=1,2 

RL3D( I ) =£QPT( BXD( I ) ) 

RLB( I ) =1o!VPLBD( 1 ) 

B (I, I )=PLB( ! ) 
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6 CONTTNMF 

C MULT S(ltJ)-^=B (I,J)=T(I,J) 

CALL GMPPD (S,B ,t,2,2,2) 

C TRANSPOSF T=TT 

CALL GMT da (T,TT,2,2) 

C MULT A’r'T = AP 

CALL GMPRD ( A , T , A P , 2 , 2 , 2 ) 
C MULT TT="AD = ALPHA 



call GMPRD (TT , AP, ALPHA,2 ,2 ,2 ) 
C FIND E VAL AND E VEC PF ALPHA 



NC = 1 

DO 21D J=l,2 
DO 21 1 1=1,2 
AX(NC)=ALPHA( I, J) 

NC=NC+1 

IF ( loGEoJ ) GO TO 212 
211 CONTINUE 
210 CONTINUE 

CALL EIGPN (AX , S , 2 , ) 
LEV=SQRT( AX(3) ) 

DO 7 1=1,2 
0 (I ) = S ( 1 , 2 ) 

7 CONTI NUF 

C MULT T7D TO GET VALUES np C 
CALL GMPRD (T,0,C,?,2,1 ) 

C EVALUATE Y AND Z 
Z(1 )=C 

DO 12 1=2,1' 

z( n=z( i-i )+^o 1 

12 CONTINUE 

DO 1? T =1 ,r 1 

Y( I )=C( 1 )7PH( 1 , I )+C(2)*PH(2, I ) 

13 CONTINUE 



WRITE (6.1 ) LEV.C r- ) ,C(2) , ( Z( I ) ,Y( I ) , 1 = 1 , 1' 1 

1C FORMAT ( •!• ////////, ’‘'' ,16X, ’ iHi- SMALLEST -IG 
lF15o6// , *f ’ , 1?X, 'THF VALUES OF Cl AND C2'//,* 
2 'Cl = ’ ,F1 5o6// , , 1 5X, »C2=' ,F15,6/ / , , 15X, 

3'THE VALUES OF Y(X) AT ONF-TENTH INTERVAL'//, 
,18X , • X* , ISX, ' Y • // ( • ’ , 2F15o 6/ / ) ) 

WRIT?- (6,9<^8) 

998 FORMAT (*1‘,2X) 

STOP 

END 
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COMPUTER DPCGRA^•i 2o 



RAYLEIGH-RITZ METHOD APPLIED TO EQN (lol) 
(1 + X)>’^Y« » (X)+V« (X)+(W**2)"Y(X)=r 



DIMENSION A(3,?) ,B(3, 3) ,AX(6> ,BX(a) , R ( 3 , 3 ) , R L B ( 3 ) , 
1T(3,3),ALPHA(3,3),AP(A,4) ,TT(3,^),S(3,3),D(^),C(3). 
2BXDC^ ) , RLBD(? ) , AO( 4 , 4 ) . BP ( ^ ^ , CM ( 4, 3 ) , C^T ( 3 t G ) » BCM( 
34,3) , ACM(4,3) ,P (ifl ) , PH(4 ,1C1 ) , PHP(4, 1 ) t FA(ir 1 ) 

4, FBI ion , Y( ! i ) , z ( in ) 

REAL LTV 
X=^' 

no 1 1 = 1,1 ■> 

C DEFINITION OF ELEMENTS OF MATRIX B 
P{ ! ) = 1. + X 
PHI 1, 1 ) = X 
PHI2 , I ) =X':'*2 

PHI3, I ) = ( I X-n /3o + ABS IX-lo /3o ) ) / 2o )=!^*2 
°HI4. I )= I I X-2o /? c + ABS I X-2o /3o ) ) /2o )**2 
C DEFINITICiN OP ELEMENTS OF matrix A 
PHPI i , 1 ) =1 
PHPI2 , 1 )=2 ^-X 

PHPI? , I ) =X-lo /-^o +ABS IX-lc /3o ) 

PHPI 4, I )=X-2o/3o+ABSI X-2o/3o ) 

92 X = X+[ o’ 1 

1 CONTINUE 
DX=Oc '^■1 

C COMPUTE AD|1,K) AND BP(I,K) 

DO 2 I =I ,4 
no 3 K=l,4 

PATNT = r ^ 

FBIMT:r' „ 

FA I 1 ) = P ( n '"PHP II , 1 ) *PHP I K , I ) 

F8I1 )=PHI 1,1 )*PHIK,i) 

DO 4 J = 2.1* 1 

FAIJ) = PIJ) 'PHPI I , J )*ohp(K, J) 

FBIJ)=PH(I , Jl'-pPIfK, J) 

FAINir = FA!NT + ( fa I J ) + F (J- 1 ) ) * DX / 2 
FBINT=FBINT+IF3IJ)+FBIJ-1) )^rX/2 
4 CONTINUE 

APII,K)=FAINT 
BPI I,K)=PBINT 
3 CONTINUE 

2 CONTINUE 
C INITIALIZE 

DO 51 KA = 1 ,4 
00 52 K3=l,3 
CMIKA , KB )=!•. ' 

52 CONTINUE 
51 CONTINUE 

C EVALUATION OF THE TRANSFORMATION MATRIX CM 
DO 53 KA=i ,3 
CMIKA,KA ) = lo' 

53 CONTINUE 

CM I 4, 1 ) =-9. 

CMI4,2 )="9 

CM (4,3) 

CALL G;'TRA(CM,CMT, 4,3) 

C MAKE TRANSFORMA"IOiN CM I TR ANSPUS ED) -'B-CM 
CALL GMPRD(P’^ CM,BCM,4,4,3) 

CALL GMPRDICMT, BCM,P,3,4,3i 
C MAKE TRANSFORMA’’ ION CM ( TR AMS POS E D ) * A'-i' CM 
r, ALL GMPRD I A R , C AC M ,4 , 4 , 3 ) 

C ALL G'nRD ( C MI , ACM , A , 3 , 4 , 3 ) 

NC = 'S 

C INITIALIZE BX ECR THE INPy^ TO EIGEN ROUTINE 
DO 2' 1 J = 1 ,3 
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COMPUTER PPOGRAM 3o 



STODOLA’S METHOD APPLIED TO EQN(lol) 
d+X)4^Y' » (X )+Y ' (X ) + ( W>!'*2 )*Y( X)=r 



C YB IS THE r-UESSED VALUE OF V AND YBN IS CO‘^PUTED VALUE OFY 
DIMENSION YEIiri ) , SS( 1 1 ) ,PI ( ICl ) ,YBN( I'd ) ,ER(1C1 ) , 
IXdCI ) 

REAL KE 

1 CONTINUE 
S = ^o 

PT( 1 )=r o 
NI = 1 
OX='^. 

BFROG=’ o 
DO 15 ] = ltKl 
X( I ) = o'‘ 1*( I-l ) 

15 CONTINUE 

DO 2 1 = 1,1" 1 
SS( I) =" 

YEU I )=-io*Xf n>d I -X( I ) ) 

C NORMALIZE YB 

IF (ABS ( YE{ I n oLEc BFROG) GO TO 2 
BFR0G=ABS(YP( ! ) ) 

2 continue 

DO 16 1=1,1"! 

YB( I I=YB( I / /PFROG 

16 CnN_i^TNi,U- 

3 CONdNUP 

DO <T 1 = 1, ICf 
SO =S 

S=S+( YP ( I } +VB ( ! + •' ) ) *DX,/R^ 

SS( I + l ) =SS ( T ) + (SO/ { • +X ( I ) )+S/( 1+X{ I + l ) ) )’^DX/2o 
PI ( I + l ) = ( 1 / n +X ( I ) ) +1-, / d+X ( I+l ) ) )=^-DX/2o +PI ( I ) 

^ CONTTN'lF 

C NOW GFT NORM OF THE NEW YB 
PFROC= o 
0051=1,11 

Y RN ( I ) = o I ( 15 ''=S F n ^ 1 ) / P T ( If. 1 ) -S S (I ) 

IF (ABS(YBN(I) )oLE,. BFROG) GO TO 5 
BFROG=ABS( Yb^'( T ) ) 

5 CONTINUE 

OMB = SORTd /BFROG) 

ERR=''o 

DO 6 1=1,1 : 

YBN( I ) =':'BN( T 5 /BF ROG 
ER( I ) = VR( n-YBN ( I ) 

IF (ABS ( ER( I ) ) o LT„ ERR 5 GO TO 6 
ERR=AB5( ER( I ) ) 

6 CONTINUE 
KE=f o 

PF = ?''o 

7 CONTINUE 

C CHECK TO SFF HOW STODOLA COMPARES WITH RAYLEIGH QUOTIENT 
DO R 1 = 1,1. ' 

Kc=KF+(YBN( T ) v:‘-?+YBM( I + 1 ) ** 2 ) - OX/ 2 

PE=PE + ( 1 + X( I + l ) ) i-M ( YBN ( I + l ) -''BN( I ) 5/DX )Y^’-2>^DX 

8 CONTINUE 
0MB2=S0RT ( PF/KE ) 

IF (MIoLTo?) go to 12 

P V;oiT[-(6,ir) NI , 0MB , err , PF , KE ,0MP2 , ( X( T ) , YRM( I ) , 

1 1=1 ,1^1 , ID) 

I*"’ FORMAT ( d'///////, ' ,15X,’NI = ' ,12, dO-- GA =',F12->6, 
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l/15Xt’THF MAXIVUr* DIFFERENCE EETWEFN YB AND YBN IS' 
2F12o6, / 13X, ' PE , KF, AND ROOT (PE/KE) AR E ' » 3F 12o 6 , / 
3/2''X,' X' ,1' X, » YBN' ,// , (12X,2F12o6/ / ) ) 

12 CONTINUE 

DO 11 1 = 1, 1^ 1 
YB( I) =YBN( I ) 

11 CONTINUE 
NI=NI+I 

IF (NIr LTo23) GO TO 27' 

GO TO 999 

270 IF ( ERPo GTo ( oO:"'! ) ) GO TO 3 
GO TO 999 
999 WRITE (6,998) 

998 FORMAT ( ' 1 ' ,2X ) 

STOP 

END 
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COMPUTER PROGRAM 4 



NEVaON'S METHOD APPLIED TO EQN (lol) 
(l + X)vY'* (X)+Y‘ (X ) + (W’^*2 )«Y(X )=0 



1 



C GET 
C USE 

ICO 



20C 



300 



DIMENSION X( 101 ) ,Y( 101) , YP( 101) , Y0( 101) , YPO( 101 ) 

DO 1 1 = 1 ,1'. 1 

Xn )=o I-l) 

CONT INUE 
Y( 1 )=0o 
YP( l)=lo 
YFN=YP( 1 ) 

DX=0oGl 
NI = 1 
YO ( 1 ) =C 
YP0( 1 ) =0 
YP0N=YP0( 1 ) 

MODE AND OMEGA BY NEWTONS METHODo FIRST GUESS OM = 

RUNGA KUTTA 2ND ORDER INTEGRATION ROUTINE 

OM=2oO 

DO 2 1=1,10.. 

J=-l 

XN=X( I+l ) 

YN=Y( I ) 

YON=YO( I ) 

YFN=YP( I )-( OM-'!''-!.2 )* ( Y( I ) / ( 1+X( I ) ) +YN/( l+XN ) )-*DX/2 
l-( YP (I ) / (1 + X( I ) ) +Y°N/ ( l+XN) )-DX/2 



YN=Y( I ) + (YP( I )+YPM)^^DX/2 

YPON-YPO ( ! ) - ( Y PO ( I ) / ( 1 +X ( I ) ) +YPnN / ( 1 +XN > ) *DX/2- ( OM 



'01 



NT 



jlV, 



-2"UM* ( Y ( X ) / 



X 



2+YN/( l + XN) )" DX/2 
YON=YO( I )+(YPO( I )+YP0N)-DX/2 
IF (JoGToO) GO TO 4..;:.' 



2o 



! .1 ) 



J=1 

GO TO 20'^ 

400 YP( I+l )=YPN 
Y(I+1)=YN 
Y0( I+l ) =YON 
YPO( I+l ) =YPON 
2 CONTINUE 

IF (NIoLToB) GO TO 11 

WRITE (6,lv ) NI ,OM, (X( I ) ,Y( I ) ,1 = 1 ,1C1 ,1( ) 

10 FORMAT (' 1 ',////////, *0 ' ,15X , 'THIS IS PATH NOo',12, 
I'o OMEGA= ' ,F12o5,//,13X,’X* ,11X,'Y' ,//(_.X,2F12o6//) 

11 CONTINUE 



NI=NI+1 



C IF ADMISSABLE QUIT 

IF (Y(lvn)oLTo (OoOCOl) ) GO TO 
C CORRECT OMEGA 

OM=OM-Y( 101 ) /Y0( 101 ) 

C IF NIoGTo25 QUIT 

IF (NIoLTo25) GO TO 100 
5G0 CONTINUE 
STOP 
END 
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COMPUTER PROGRAM 5o 
RAYLEIGH-RITZ METHOD APPLIED TO 
BESSEL'S EQUATION OF ORDER ZERO 



DIMENSION A(3, 3) ,6( 3,3) ,AX(6) ,BX(6) , R ( 3 , 3 ) , RLB { 3 ) 

1T(3,3),ALPHA(3,3),AP(4,^),TT(3,3),S(3,3),D{3),C(3), 
2BXD(3) ,PLBD(3) ,AQ(4,4) ,BP(4,A) ,CM{4,3) ,CMT(3,4), 
3BCM(4, 3 ) ,ACM(4,3) , P ( 101 ) ,PH{ 4, ICl ) ,PHP(4, 101 ) ,FA( 101 ) 
4FB(101) ,EM(101 ),Y{1C1) tZdCl) 

REAL LEV 
X=0o 

DO 1 1=1,101 
P( I ) = X 
EM{ I )=X 

C DEFINITION OF ELEMENTS OF MATRIX B 
PH( 1, I )=lo 
PH(2, I )=X**2 

PH(3,I ) = ( ( X-oB + ABS ( X-o 8 ) )/2o )=!^-2 
PH(4, I )=(( X-lo6+A8S(X-l«6) )/2o )=^-'*2 
C DEFINITION OF ELEMENTS OF MATRIX A 
PHP{ 1, I )=2o-X 
PHP( 2, I )=X-o8+ABS( X-o 8 ) 

PHP{3. I )=X-lo6+ABS(X-lo6) 

92 X=X+Co' 240 5 

1 CONTINUE 
DX=(.o02AO5 

C COMPUTE AP(I,K) AND BP(I,K) 

DO 2 1=1,4 
DO 3 K=l,4 
FB INT=' o 

FB{1) = EM(1 )'^PH( I,i)*DH(K,l) 

DO 4 J = 2.1.:'] 

FBI J ) =EM( J )*PH ^ T . j ) ;.= PH I K, i ' 

FBINT=FBTNT+(F8( J ) +FB\ J-i ) )*DX/2 
4 CONTINUE 

BP( I,K )=FBINT 
3 CONTINUE 

2 CONTINUE 

DO 600 1=1,3 
DO 601 K=l,3 
FAINT='' o 

FA(1 ) = P(1)'-PHP{ I,1)--!=PHP(K,1) 

DO 602 J=2,iCi 

FA( J) = P( J)--PHP( I ,J)*PHP(K, J) 

FAINT=FAINT+(FA( J)+FA( J-1) )*DX/2 
602 CONTINUE 

API I,K )=FAIKT 
AII,K)=API I,K) 

601 CONTINUE 
600 CONTINUE 
C INITIALIZE CM 

DO 51 KA=1 ,4 
DO 52 KB=1,3 
CMIKA,KB)=OoC 

52 CONTINUE 
51 CONTINUE 

C EVALUATION OF THE TRANSFORMATION MATRIX CM 
DO 53 KA=1,3 
CMIKA + 1 ,KA ) = lo0 

53 CONTINUE 

CM! 1, 1 ) =-5o784 
CMIl,2)=-2o576 
CM! 1 ,3) =-Oo 648 
CALL GMTRA (CM,CMT, 4, 3) 

C MAKE TRANSFORMATION CM I TRANS POS E D ) * 6- CM 
CALL GMPRD(BP,CM,8CM,4,4,5) 

CALL GMPRDICf'T, BCM,6,3,4,3) 

NC = 1 
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C INITIALIZE BX FOR THE INPUT TO EIGEN ROUTINE 
DO 2P1 J=1 ,3 
DO 202 1=1,3 
BX(NC)=B(I , J) 

NC=NC+1 

IF { loGEoJ ) GO TO 201 
202 CONTINUE 
2C1 CONTINUE 

CALL EIGEN (BX,R,3,0) 

IF(BX(6)ol.ro' o ) GO TO 80 6 
C RLB(I) IS 1/SQRT OF Eo VALo OF B 
KK=1 

DO 8 11=1,3 
BXD( I I )=BX( KK) 

KK=KK+I I+l 

8 CONTINUE 

DO 1C 11=1,3 
DO 9 JJ=1,3 
B (II,JJ)=Co 

9 CONTINUE 
10 CONTINUE 

DO 6 I =1,3 

RLBD( I )=SQRT(BXD( I ) ) 

RLB( I )=loC/RLBD( I > 

B (1,1) =RLB( I) 

6 CONTINUE 

C MULTIPLY R(I,J)-8 (I,J)=T(I,J> 

CALL GHPRD(R,B ,T,3,3,3> 

C TRANSPOSE T=TT 

CALL GMTRA( T,TT ,3,3 ) 

C MULTIPLY A'^=AQ 

CALL GMPRD( A,T, A0,3,3,3) 

C MULTIPLY TT*AQ=ALPHA 

CALL GMPRD( TT,AO,ALPHA,3,3,3) 

C FIND EoVALo AND EoVECo OF ALPHA 

Mr Li ’ - - • 

— A. 

DO 210 J=l,3 
DO 211 1=1,3 
AX(NC)=ALPHA(I , J) 

NC=NC+1 

IF( loGEo J ) GO TO 210 
211 CONTINUE 
210 CONTINUE 

CALL EIGEN (AX ,S,3,0) 

C LEV IS THE LOWEST EIGIN VALUE 
LEV=SQRT( AX(6) ) 

DO 7 1=1 ,3 
D( I >=S( I ,3) 

7 CONTINUE 

C MULTIPLY T>)=D TO GET VALUES OF C 
CALL GHPRD(T,D,C,3,3,1) 

C0NE=-5o76A^C( 1 ) -2 . 3 76^C ( 2 ) -( o 6AG^C ( 3 ) 

z( n=0ov 

DO 64 1=2,17 1 

Z( I )=0o024 '5^( I-l ) 

64 CONTINUE 

DO 65 1=1, U I 

Y( I )=C0NE^:<PH(1 , I )+C (1 ) *PH(2, I ) +C ( 2 ) *PH( 3 , I ) +C ( 3 ) * 
1PH(4,I ) 

65 CONTINUE 

WRITE ( 6, ICQ )LEV,CONE,C ( 1) ,C(2),C(3),(Z(I),Y(I), 
11=1,101 , 1C) 

109 FORMAT (' 1 '////////,' C ', 1 5 X ,' THE SMALLEST EIGENVALUE', 
lF15o6//, 'O' ,15X, 'VALUES OF C1,C2,C3,AND C4 1 5X 
2, *C1= • ,F15o6//, ,15X, 'C2=' ,F15o6//,'0' ,15X, 'C3 = ' , 
3F15o6//,'0',15X, «C4= ' , F 1 5 , 6// , • C • , 1 5 X , 

4'VALUES OF Y(X) AT ONE-TENTH INTERVAL'//, 

5 'O' , 19X, 'X' , 14X, ' Y' //( ' ' ,2F15c6// ) ) 

806 STOP 
END 
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COMPUTER PROGRAM h 



STODOLA'S method APPLIED TO 
BESSEL'S EQUATION OF ORDER ZERO 



C YB IS THE GUESSED VALUE OF Y AND YBN IS CO'^PUTEO VALUE OFY 
DIMENSION YCn , SS( l‘n ) ^ PI ( I'U ) ,vbN(] " 1 ) ,ER( 1' 1 ) , 

IX ( 101 ) 

REAL KE 

1 CONTINUE 
S=^o 
MI=1 

DX=Cc CPA'^5 
BFROG=' o 
DO 15 I=1,1M 
X( I )=r, 024 5=;:( l-l ) 

15 CONTINUE 

DO 2 I=1,1E‘' 

SS( T ) =” 

YR( I }= ( 1-X( I )/2.4'-E 
C NORMALIZE YB 

IP ( ABS( v8 ( I ) K lEo EFROG) GO TO 2 
BFROG= ABS(Ve{ I ) ) 

2 CONTINUE 

DO 16 1 = 1,1- 1 
YB( I )=YR( I )/BFRDG 

16 CONTINUE 

3 CUNTINUF 

SSI • != ^ 

SS(?)=YB(2)=-DX/2 
S = * o 

S = S+(YB( 1 )=:'«' (1 )+X(? )*YB(2) )#DX/2 
DO 4 I=2,K- 
SO =S 

s=s+('^'p ( n-x( I ) +Y( 1+1 ) -!=vb( 1 + 1 ) )*nx /2 
SS( I + l ) =SS ( I) + ( so/ X( I ) +S/X( I + l ) )->DX/2 
^ CONTINUE 

C NOW GET NORM OF THF NEW YB 
BFROG = ' o 
DO 5 1 = 1 r 1' 1 
VBN( I 5 =SS( i ' ] )-SS ( I ) 

IF (A8S(YBN( I ) )oLE,BFROG) GO TO 5 
BFROG=ABS(YPN( I ) ) 

5 continue 

0M5 = SQPT( 1 . /BFPOG) 

ERR=' o 

DO 6 1=1,1^ 1 

Y6N( 1 ) = Y8N( I ) /RFP.OG 

ER( I )=Y6( I )-YBM( I ) 

IF (ABS(ER( I) loLToEPR) 'GO TO 6 
ERR=A3S f ER ( I ) ) 

6 CONTINUE 
KE = '' o 

P E = '' o 

7 CONTINUE 

C CHECK TO SEE HOW STODOLA COMPARES WITH RAYLEIGH QUOTIENT 
DO 8 1=1,10' 

KE = KE + ( X( I YPN ( I )* ■'2 + X ( I + l )*Y3N( 1 + 1 )'^-^2 ) =^=DX/2 
PP = PF +X( I ) = ( ( VBNH + 1) - YBN( I ) ) /DX)'J'"2''0X 

8 CONTINUE 
0MB2=S0RT(PF/KP) 

IF (NIoLToT) GO TP 1? 

o WRITF(6,10) NI ,0WB. pRR, PP, KE,0MB2, ( X( I ) , ^’ '(I), 

11=1,101 ,if') 
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FORMAT (' 1 '///////,' •, = S I? O OMEGA= F12o 6, 

1/15X,»THE 'MAXIMUM DIFFERENCE BETWFpN Y8 AND YBN IS' 
2F12o6, /] 5X, ' RE, KEt AND ROOT (PE/KE) AR E ' , 3^ 12„ 6, / 
3/2':X,»X' ,10X, 'YBN' ,//, (l2X,2F12o6//) ) 

12 CONTINUE 

DO 11 I =1,1' 1 
YB( I )=Y6N( I ) 

11 continue 

NI=NI+1 

IF (MIuLTo25) GO TO 27'R 
GO TO 999 

270 IF (FRRoGTo (orTUl ) ) GO TO 3 
GO TO 090 
999 write (6,99S) 

998 FORMAT ( » 1’ ,2X) 

STOP 

END 
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COMPUTER PROGRAM 7 



NEWTON'S METHOD APPLIED TO 
BESSEL'S EQUATION QR ORDER ZERO 



DIMENSION Xn/'l ) ,v{ 1* 1 ) , YP( 1' 1 ) ,YO( 1^1 ) , YPO( l''I ) 
1,Y01(1^ 1) ,Y02( IPl) 

DO 1 1=1,11 

X( I )=0 024REv-( l-l ) 

1 CONTINUE 
Y ( 1 ) = 1 ' 

YPM ) 

YPN=YP( 1 ) 

DX=f\,D?4 5 
NI=1 

YG( 1 )=V, 

VPOd )=''o 
YPnN=YPO( 1 1 

C GET MODE AMD OMEGA BY NEWTONS METHODo FIRST GUESS CM = 2or 
C USE RUNGA KUTTA 2ND ORDER INTEGRATION ROUTINE 

Of‘ = lo4 

ICO XN=X(2) 

YM=Y( 1 > 

YQN = YO( 1 ) 

YPN=YP( 1 )-(OM^v2)*(Y(l )+YN)*DX/2o-(-o?-(OM**2)’’=Y( 1) 

1 -o5^^0M=!^^2)'’ DX/2o 
1 = 1 

YN=Y ( 1 ) +( YPd ) +YPN ) ^nX/2o 

YPN=YP( 1 )-( OM*=!=2 )* ( Yd )+YN) *DX/2 - - (-o 5* ( OM*^'=? ) *Y ( ] ) 
1+YPN/XN ) >'DX/?o 
YN=Y (1 ) + ( YPd ) +VPN )'dX/2 

YP0N=YP0( 1 »-( -OM>-Yn I-OM) =t^0X/2o -(OM>^'’‘2 \YOd )+YCN)--nX 
1 /2o -2 o ( Y< I j+V'N i •“!>'//- 

YON=YOd ) + (vr 0(1 ) + v pqm ) . nX/2 ~ 

YPON=YPO( 1 ) - ( - 0M=‘ Y( 1) + YP0N/Xfd*DX/2o- (0MY-*2 ) V- (YOd ) 

1 +YON )*DX/2o *0M-'= (Yd )+YM ) ^"DX/2-, 

YON=YO( 1 ) + (YPOd ) +YP0N)v-DX/2 

YP( 2)=YPM 

Y(?)=YN 

Y0( 2) =YON 

YPO( 2 ) = YPOI\’ 

no 2 1 = 2 , U 

J=-l 

YN=Y( I ) 

YON=VO( I ) 

2S-D YPN = YP( I ) - ? ) *( Y( T ) + YN) -^DX/2o -(vp ( I ) / X ( I ) +YPN/XM)=!-- 

lDX/2 

YN=Y( I) + ( YP( I ) + YPN)*nx/:>. 

YPON=YPO( I ) - ( YPOd ) / v ( n +YPOr /XN)=rDX/2o -( 0M*'*2 )=^- ( Y0( I) 
1+YON)=7DY/2 o-2o-OM’’= (Y( I )+VN ):;^ox/2" 

Y0N=V0( I ) + (YPO( I )+vpon)>-ox/2c 
IP doGTo5) GO TO 54 
54 CONTINUE 
300 IE (JoGTc^') GO TO 
XN=X( I+l ) 

J = l 

GO TO ?r 

40P YP(I+1)=VPN 
Y( 1+1 ) =YN 
Y0( I+l ) =Y0N 
vpo( 1+1 )=YPGN 

2 CONTINUE 

IF (NIoLTo-'O 3 0 TO 11 

WRITE (6,K) NT T CM, (X ( I ) , Y ( I ) , 1 = 1 , ir 1 , 1 '' ) 

1'' FORMAT ('I'-////////,’ ',15X,’THIS IS PATH NCo',12, 

Id OMEGA= ' , F12. 5, //, 1 3X , • X' , IIX, » V. , //( ’ ' x,2F12o 6/ /) ) 
1] CONTINUE 
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NI=NI+i 

C IF ADMISSIBLE DUI 

IF (ABS(V(l']))oLTo (<^oOCCOI ) ) GO TO 15 
C CORRECT OMFGA 

OM=OM-Y( 1C3 ) /YO( l'* l ) 

C IF MIoGTc25 QUI" 

IF (NIo LTo25) GO TO ICO 
15 CONTINUE 
STOP 
END 
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